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Abstract 

We present a theoretical and numerical framework to compute bifurcations of 
equilibria and stability of slender clastic rods. The 3D kinematics of the rod 
is treated in a geometrically exact way by parameterizing the position of the 
centerline and making use of quaternions to represent the orientation of the ma- 
terial frame. The equilibrium equations and the stability of their solutions are 
derived from the mechanical energy which takes into account the contributions 
due to internal moments (bending and twist), external forces and torques. Our 
use of quaternions allows for the equilibrium equations to be written in a simple 
quadratic form and solved efficiently with an asymptotic numerical continua- 
tion method. This finite element perturbation method gives interactive access 
to semi-analytical equilibrium branches, in contrast with the individual solution 
points obtained from classical minimization or predictor-corrector techniques. 
By way of example, we apply our numerics to address the specific problem of 
a naturally curved rod under extreme twisting and perform a detailed com- 
parison against our own precision model experiments of this system. Excellent 
quantitative agreement is found between experiments and simulations for the 
underlying 3D buckling instabilities and the characterization of the resulting 
complex configurations. 
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1. Introduction 

Filaments, rods and cables are encountered over a wide range of length- 
scales, both in nature and technology, providing outstanding kinematic free- 
dom for practical applications. Given their slender geometry, they can undergo 
large deformations and exhibit complex mechanical behavior including buckling, 
snap-through and localization. A predictive understanding of the mechanics of 
thin rods has therefore long motivated a large body of theoretical and computa- 
tional work, from Euler's elastica in 1744 [1] and Kirchhoff's kinetic analogy in 
1859 [2] to the burgeoning of numerical approaches such as finite clement-based 
methods in the late 20 th century [3], and the more recent algorithms based on 
discrete differential geometry [4] . Today, these advances in modeling of the me- 
chanics of slender elastic rods are helping to tackle many cutting-edge research 
problems. To name just a few, these range from the supercoiling of DNA [5], 
self-assembly of rod-coil block copolymers [6] , design of nano-electromechanical 
resonators [7, 8], development of stretchable electronics [9], computed animation 
of hairs [10] and coiled tubing operations in the oil-gas industries [11]. 

An ongoing challenge in addressing these various problems involves the ca- 
pability to numerically capture their intrinsic geometric nonlinearities in a pre- 
dictive and efficient way. These nonlinear kinematic effects arise from the large 
displacements and rotations of the slender structure, even if its material prop- 
erties remain linear throughout the process [12]. As a slender elastic rod is 
progressively deformed, the nonlinearities of the underlying equilibrium equa- 
tions become increasingly stronger leading to higher densities in the landscape 
of possible solutions for a particular set of control parameters. When multiple 
stable states coexist, classic step-by-step algorithms such as Newton-Raphson 
methods [13] or standard minimization techniques [14] are often inappropriate 
since, depending on the initial guess, they may not converge towards the de- 
sired solution, or any solution. Addressing these computational difficulties calls 
for alternative numerical techniques, such as well-known continuation methods 
[15, 16]. Continuation techniques are based on coupling nonlinear algorithms 
(e.g. predictor-corrector [15] or perturbation methods [16]) with an arc-length 
description to numerically follow the fixed points of the equilibrium equations 
as a function of a control parameter, that is often a mechanical or geometrical 
variable of the problem. With the goal of determining the complete bifurca- 
tion diagram of the system, these methods enable the computation of all of the 
equilibrium solution branches, as well as their local stability. 

Two main approaches can be distinguished for continuing the numerical 
solutions of geometrically nonlinear problems. The first includes predictor- 
corrector methods whose principle is to follow the nonlinear solution branch 
in a stepwise manner, via a succession of linearizations and iterations to achieve 
equilibrium [13]. These methods are now widely used, particularly for the nu- 
merical investigation of solutions of conservative dynamical systems, with the 
free path-following mathematical software AUTO being an archetypal exam- 
ple [17]. Quasi-static deformations of slender elastic rods have been intensively 
studied using this software [18, 19], mostly due to the analogy between the rod's 
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equilibrium equations with the spinning top's dynamic equations [20]. Although 
popular and widely used, the main difficulty with these algorithms involves the 
determination of an appropriate arc-length step size, which is fixed a priori 
by the user, but may be intricately dependent on the system's nonlinearities 
along the bifurcation diagram. A smaller step size will favor the computation of 
the highly nonlinear part of the equilibrium branch, such as bifurcation points, 
but may also impractically increase the overall computational time. On the 
other hand, a larger step size may significantly compromise the accuracy and 
resolution of the results. 

The second class of continuation algorithms, which have received less at- 
tention, is a perturbation technique called the Asymptotic Numerical Method 
(ANM), which was first introduced in the early 1990's [21, 22]. The underly- 
ing principle is to follow a nonlinear solution branch by applying the ANM in 
a stepwise manner and represent the solution by a succession of local polyno- 
mial approximations. This numerical method is a combination of asymptotic 
expansions and finite element calculations which allows for the determination of 
an extended portion of a nonlinear branch at each step, by inverting a unique 
stiffness matrix. This continuation technique is significantly more efficient than 
classical predictor-corrector schemes. Moreover, by taking advantage of the an- 
alytical representation of the branch within each step, it is highly robust and 
can be made fully automatic. Unlike incremental-iterative techniques, the arc- 
length step size in ANM is adaptative since it is determined a posteriori by 
the algorithm. As a result, bifurcation diagrams can be naturally computed in 
an optimal number of iterations. The method has been successively applied to 
nonlinear elastic structures such as beams, plates and shells but the geometri- 
cal formulations were limited to the early post-buckling regime and to date, no 
stability analyses were performed with ANM [23, 24, 25]. 

In this paper, we develop a novel implementation of the semi-analytical ANM 
algorithm to follow the equilibrium branches and local stability of slender elastic 
rods with a geometrically-exact 3D kinematics. In Section 2, we first describe the 
3D kinematics where the rod is represented by the position of its centerline and 
a set of unit quaternions to represent the orientation of the material frame. In 
Section 3, we then derive the closed form of the rod's nonlinear equilibrium equa- 
tions, by minimizing its geometrically-constrained mechanical energy including 
internal bending and twisting energy, as well as the work of external forces and 
moments. Expressing the flexural and torsional internal moments in a quater- 
nion basis yields differential equilibrium equations that are simply quadratic in 
terms of the unknowns. In Section 4, we proceed by presenting the numerical 
method developed to compute the equilibrium solutions. Using a finite element 
approach, the discretized system of equilibrium equations can be solved with 
the ANM algorithm, which is particularly efficient for computing our algebraic 
quadratic form. The local stability of the computed equilibrium branches is 
assessed by a second order condition on the constrained energy. Finally, we 
describe how to implement this numerical method in the open source software 
MANlab; a user-friendly, interactive and Matlab-based path-following and bi- 
furcation analysis program [26, 27]. In Section 5, we develop our own precision 
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model experiment for the fundamental problem of the writhing of a clamped 
elastic rod [18, 28, 29, 30] and challenge our numerical model against experi- 
mental results. The simulations are robust, computationally time-efficient and 
exhibit excellent quantitative agreements with our experiments, demonstrating 
the predictive power of our framework. 

2. Kinematics 

In this section, we present the formulation for the geometry and 3D kine- 
matics of the slender elastic rod that we will use in our study. Assuming no 
shear strains and inextensibility, the mechanical deformations are represented 
by the rate of change of the orientation along the rod, characterized by a set of 
geometrically constrained unit quaternions. 

2.1. Cosserat theory of elastic rods 

An elastic rod is a slender elastic body which has a length along one spatial 
direction that is much larger than its dimensions in the two other perpendicular 
directions, that define the cross section [Fig. 1(a)]. We denote the typical size 




Figure 1: Kinematics of the Cosserat rod in the global cartesian frame (x,y,z). (a) The 
configuration of the rod is defined by its centerline r(s). The orientation of each mass point of 
the rod is represented by an orthonormal basis (d\(s) d,2(s) ££3(3)), called the directors, where 
£(3(5) is constrained to be tangent to r(s). (b) The three local modes of deformation of the 
elastic rod, associated with the change of (bl) material curvature rei related to the direction 
d\ of the cross-section, (b2) material curvature ft 2 related to d.2, and (b3) twist. 
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of the cross-section by h and the other length scale by L. At large scales, the 
rod can be regarded as an adapted material curve: its centerline. If s denotes 
the curvilinear coordinate along the centerline of the undeformed rod, we can 
represent this line by a position vector function (with respect to some fixed 
origin) of the material point originally at s in the reference configuration, 

r(s) = [r x {s) r y (s) r z (s)} T = [x(s) y(s) z(s)f . (1) 

We consider unstretchable rods whose centerline remains inextensible upon de- 
formation. As explained in detail in [12], this assumption is physically justified 
for a wide range of loading conditions, provided that the aspect ratio of the 
rod, h/L, is small. Under this assumption, the variable s is also the curvilinear 
coordinate along the centerline in the actual configuration. The configuration 
of the rod is not only characterized by the path of its centerline but also by how 
much it twists around this line. We consider this twist by introducing the mate- 
rial frame (di(s) d 2 {s) d 3 (s)) in the deformed configuration. At each particular 
location s, we associate an orthonormal basis dk, (k = 1,2,3) attached to the 
centerline. The centerline, together with this set of material frames, form what is 
called a Cosserat curve [31]. We choose the orientation of these material frames 
in a way such that the directors d\ and d 2 lie in the plane of the cross-section, 
while the third director d 3 is always parallel to the tangent of the curve [see 
Fig. 1(a)]. Considering the case of small strains, the triad (di(s) d2{s) d 3 (s)) 
remains approximately orthonormal upon deformation. This is known as the 
Euler-Bcrnoulli kinematical hypothesis (assumption of no shear deformations). 

Before we are able to establish the constitutive relation, we have to quan- 
tify the rate of change of position and orientation along the rod's centerline. 
The rate of change in the position of the centerline is a strain vector v(s) = 
[v i(s) v 2 (s) V3(s)] that vanishes since shearing in both transverse directions and 
stretching are neglected. Therefore, the strains arise from the orientational rate 
of change of the cross-sections alone, which we now express using the framework 
of differential geometry of curves in 3D space [12]. The previous condition of 
orthonormality (Euler-Bernoulli assumption) yields the relations, 

<(s).di(s) = and d'^.d^s) = -t^(s).d;(s), (2) 

for all indices i and j varying from 1 to 3 (there is no implicit sum over i in the 
first equation) and where ( )' denotes differentiation with respect to s. The most 
general set of first-order linear equations conserving the orthonormal character 
of the material frame (di(s) d 2 (s) d 3 (s)) represented in Eq. (2) can be expressed 
as, 

d[(s) = K 3 (s)d 2 (s) - K 2 (s)d 3 (s), (3a) 
d 2( s ) = -K S (s)d 1 (s) + K 1 (s)d 3 (s), (3b) 
d' 3 (s) = K 2 (s)d 1 (s) - K 1 (s)d 2 (s). (3c) 
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where Ki(s), K2(s) and k 3 (s) are scalar functions interpreted below. These 
equations describe the rigid-body rotation of the frame. Using the notation 
u x v for the cross product of two vectors, Eqs. (3) can be rewritten as, 

di(s) = n(s) x di(s), d' 2 {s) = tl(s) x d 2 {s), d' 3 (s) = n(s)xd 3 (s), (4) 

where we have introduced the Darboux vector J~2(s), 

U(s) = K 1 (s)d 1 (s) + K 2 (s)d 2 (s) + K 3 (s)d 3 (s). (5) 

The physical interpretation of Eqs. (4) is that the material frame rotates with 
a rotation velocity, Jl(s), when following the centerline at unit speed. The 
quantities k\ and K2 in Eq. (5), called the material curvatures, illustrated in 
Fig. I(bl)-(b2), represent the extent of rotation of the material frame, with re- 
spect to the directions d\ and d 2 of the cross-section. The quantity n 3 quantifies 
the rotation of the material frame with respect to the tangent d 3l and is called 
the material twist of the rod [see Fig. I(b3)]. In order to write the material 
curvature and twist in an explicit form, the Darboux vector has to be rotated 
into the local frame. Using the condensed notation, k(s) — \k\{s) k 2 {s) k 3 (s)] T , 
and the rotation matrix of the Euclidean 3D space R(s) €K 3x3 , this rotated 
Darboux vector is, 

#c(s) = R T (s)n{s), (6) 
or, in terms of the directors dfe(s), 

= d k (s).ft(s), (7) 

since the directors dfc(s) constitute the columns of the rotation matrix R(s) — 
[di(s) d 2 (s) d 3 (s)]. 

The 3D kinematics formulation of our inextensible and unshearable elastic 
rod is not yet complete because we won't be able to derive the equilibrium 
equations directly from the material curvatures. In fact, a difficulty arises when 
trying to compute the infinitesimal work of the external forces using the variables 
Ki, k 2 and K3. A perturbation of these quantities yields a non-local perturbation 
to the centerline and attached material frame so that the work of the external 
forces cannot be written in a straightforward manner [12]. Instead, the classic 
approach is to choose as degrees of freedom the orientation of the material frame 
characterized, in this paper, by a set of quaternions. We shall now explain how 
to represent the rotation matrix R(s) or the directors, dk(s), and the strain 
rate vector, k(s), in the framework of quaternions. 

2.2. Quaternion representation 

Quaternions are a number system that extends the complex number repre- 
sentation of geometry in a plane to the three-dimensional space [32] . They were 
first described by Hamilton in 1843 [33, 34] and were extensively used in many 
physics and geometry problems before loosing prominence in the late 19 th cen- 
tury following the development of numerical analysis. Quaternions were then 



6 




Figure 2: Rotation of a rigid body using Euler's rotation theorem and a set of unit quater- 
nions. Knowing the rotation angle <f> around the unit vector 6, we can associate a rotation 
matrix R(q) and three directors di(q), d2[q) and d^{q) expressed exclusively in terms of 
quaternions according to Eqs. (10)-(13). 

revived in the late 20 th century, primarily due to their power and simplicity 
in describing spatial rotations, and have since been revived in a wide range 
of fields: applied mathematics [35], computer graphics [36, 37], optics [38, 39], 
robotics [40] and orbital mechanics [41, 42]. It is beyond the scope of this article 
to discuss a detailed evaluation of the advantages and disadvantages of using 
quaternions over other rotation parametrizations. However, we highlight that 
quaternions are a non-singular representation of rotation, unlike Euler angles 
for instance, even if they are less intuitive than direct angles. Moreover, we 
favor quaternions over trigonometric approaches because of their remarkably 
compact quadratic polynomial form. We will show that one striking outcome of 
using quaternions is that the equilibrium equations we shall derive are, at most, 
cubic in terms of the degrees of freedom. This property is at the heart of the 
numerical continuation method presented in Section 4. 

The fundamental relation of the algebra of quaternions, denoted by W, is, 

i 2 = j 2 = k 2 = ijk = -1, (8) 

where i, j, and k are the basis elements of T-L. A quaternion number in H is 
written in the form qii + (?2j + 53^ + <?4 where the imaginary part q%i + q2.j + q?,k 
is an element of the vector space 1Z 3 and the real part q^ is a scalar. Using 
the basis i, j, k, 1 of % makes it possible to write a quaternion as a set of 
quadruples, usually expressed as a vector in T? 4 , 

q = [<7i<72<?3<74] T - (9) 

Quaternions of norm one, or unit quaternions, are a particularly convenient 
mathematical notation for representing orientations of objects in three dimen- 
sions. Using Euler's rotation theorem which states that a general re-orientation 
of a rigid-body can be accomplished by a single rotation about some fixed axis, 
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one can represent a rotation by a set of quaternions, known as Euler parameters, 



q = [& s sin($/2) 6 a sin($/2) 6 2 sin($/2) cos($/2)] T , (10) 

where $ is the Euler principal angle and b = [b x b y bz] T is the unit length princi- 
pal vector such that b x +by+b^ — 1 [see Fig. 2] . Given that four Euler parameters 
are needed to define a three-dimensional rotation, a natural constraint equation 
prescribing that q is indeed a unit quaternion follows from Eq. (10), 

£ + + = 1- (11) 

The orthogonal matrix representation corresponding to a rotation by the quater- 
nion q = qii + q 2 j + qzk + q 4 with ||g|| = 1 is, 

9l - 0$ ~ ?! + <& 2(qi<2<2 - 9394) 2(gi<? 3 + 17294) 

R(Q) = 2 (<7i92 + 9394) -q\ + q 2 - ql + q% 2(g 2 93 - 9194) 

2(9193-9294) 2(q 2 93 + 9194) -9i- 92 + 93 + 94 

(12) 

Returning to the context of a thin elastic rod discussed above, its material 
frame (di(s) d 2 (s) d^{s)) remains orthonormal upon deformations and its rigid- 
body re-orientation can be expressed by the rotation matrix given in Eq. (12) 
such that, 

R(q(s)) = [di(q(a)) d 2 (q(s)) d 3 (q(s))} . (13) 

The local frame is now parametrized in terms of the curvilinear unit quaternion 
coordinates vector q(s) — [qi(s) 92(5) 93(s) 94(s)] T along the slender rod. 

We proceed by relating the strain rate vector k(s) of Eq. (6) to the Euler 
parameters q(s). Multiplying each of the three geometric relations given in 
Eqs. (3) by the relevant director d k yields expressions for the material curvatures 
and twist in terms of the directors alone, 

Ki(s) = -d 2 (s).d' 3 {s), k 2 (s) = -d 3 (s).di(s), « 3 (s) = -di(s).d' 2 (s). (14) 

To compute d^(s) in terms of quaternions, we note that d' k (s) is a function of 
q(s), which is itself a function of the curvilinear coordinate s. Upon employing 
the chain rule of partial differentiation, we obtain 

j/ 1 \ jt f t \\ dd k {q(s)) dd k dq , 

dfc(s) = d k ( q (s)) = — — — = ~q^~q~ s = J k(q(s))q (*) (15) 

for the three directions k — 1,2,3, where Jf.(q(s)) is the Jacobi matrix = 
ddk/dq. Replacing dfe(s) and d' k (s) by their respective expressions (12)-(13) 
and (15) in Eq. (14) allows us to express the material curvatures, Ki(s) and 
k 2 (s), and the twist, k 3 (s), solely in terms of the unit quaternions, 

K fc (s) = 2B k q( S )q'{s) fork=l,2,3, (16) 
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where the skew-symmetric matrices £?& read 

10 
-10 
1 
0-10 
(17) 

With the new expression of Kfc(s) hi Eq. (16), we have been able to write the 
total strains in terms of the locally perturbable variables q(s), which will be 
used in the derivation of the equations of equilibrium by variation of clastic 
energy presented below in Section 3. 

It is important to note, however, that the kinematic formulation is not yet 
complete since the four quaternions qi(s), ?2(s), qa (s) and <?4(s) are not geo- 
metrically independent. First, to represent a three-dimensional rotation with 
four coordinates, the unit quaternion assumption ||q(s)|| = 1 given in Eq. (11), 
must be verified. Secondly, whereas thus far we have treated the centerline po- 
sition r(s) and the orientations q(s) as separate entities, the positions and the 
orientations cannot be considered independently. Indeed, the material frames 
parametrized by r(s) and q(s) are coupled by the constraint that the third 
director d 3 (q(s)) is always parallel to the tangent r'(s), 

r'(s) = d 3 (q(s)), (18) 

where r'(s) is the unit tangent vector to the Cosserat curve and ||r'(s)|| = 1 
along the centerline since we assumed inextensibility. The three constraints set 
by Eq. (18) assure that the directors are adapted to the Cosserat curve [see 
Fig. 1]. 

The three-dimensional kinematics of our inextensible and unshearable rod 
(including bending and twist) is represented by Eq. (16), which links the strain 
rates to the local orientation of the material frame, together with the four geo- 
metrical constraints given in Eq. (11) and Eq. (18). For the remainder of this 
article, the three positions r x (s),r y (s),r z (s) and four quaternion coordinates 
Qi( s )y 92(s)) ?3(s), 94(5) constitute the seven degrees of freedom of our slender 
elastic rod [31, 43]. After taking into account the four constraint equations, 
only three of the DOFs are, in fact, geometrically independent. Their values 
are determined by the three-dimensional equilibrium equations, which we now 
address in the following section. 

3. Mechanical equilibrium 

Having formulated the kinematics of our system, we proceed by analyzing the 
energetics of an arbitrary configuration of the slender elastic rod. We will then 
derive the equations for equilibrium obtained under the assumption that this 
energy is stationary under small deformations for the given boundary conditions 
and geometrical constraints introduced above. We highlight the fact that the 
equilibrium equations are highly nonlinear due to geometry, rather than the 
material response. 
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3.1. Energy formulation 

For simplicity, and to avoid loss of generality, we shall adopt the framework of 
Hookean elasticity and consider linear isotropic constitutive laws. For practical 
purposes, this hypothesis is usually appropriate since, for slender elastic rods, 
the strains at the material level are typically small. Under this assumption, the 
total elastic energy of the slender elastic rod can be written as the uncoupled sum 
of bending and twisting contributions [12]. Although the reference configuration 
of the rod is assumed to be stress-free, we can readily account for rods with 
intrinsic natural curvature and twist. Doing so, the elastic energy of a rod with 
length L and a constant cross-section reads, 



where we used the previously defined rotational strain rate vector k(s) = 
[ki(s) k 2 (s) «3(s)] , and where the quantities «i(s), and «3(s) are the 

intrinsic natural curvature and twist of the rod along the directors d\, d 2 and 
££3, respectively. In this expression, E is the Young's modulus of the material 
and G = E/2(l + v) is the shear modulus of the material with Poisson's ratio 
v. The constants 1\ and I 2 are the moments of inertia along the principal di- 
rections of curvature in the plane of the cross-section d\ and d 2 and J is the 
moment of twist which, similarly to I\ and I 2 for the bending energy, depends 
only of the geometry of the cross-section. Replacing the material curvatures 
Ki(s), k 2 (s) and twist n^{s) by their expression given in Eq. (16) allows us to 
write the elastic energy £ e in a more compact form, in term of the rotational 
degrees of freedom q{s) alone, 



where E\ — E 2 — E, E 3 = G and I3 = J. 

3.2. Variation of the energy 

We now follow a variational approach for the elastic energy in Eq. (20), 
and consider an infinitesimal perturbation from an arbitrary configuration of 
the rod. The perturbed quantities are preceded by 5. Carrying out the first 
variation of Eq. (20), the corresponding variation of the energy £ e is, 



where 8q = \8q\ 5q 2 Sq^ Sq^] is the vector of the arbitrary perturbations of 
the rotational degrees of freedom q. Upon integration by parts, we transform 





(20) 



S£ e = E k I k / (2B k qq' - k k ) (2B k Sqq' + 2B k q8q') ds, (21) 
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Eq. (21) into an integral that depends on Sq alone to arrive at, 



S£ e 



J2 E kh (2B k qq' - k k ) 2B k qSq 



fc=i 



r L 3 

J ° k=i 



[(2B k qq' - K k ) AB k q' + {2B k q'q' + 2B k qq" - k' k ) 2B k q] Sqds, (22) 



where the first term stands for the variation of elastic energy over the entire 
interval and is the boundary term from the integration by parts assuming that 
the rod is parametrized from s = to s = L. Physically, this first term repre- 
sents the work done by the operator upon a change of orientation applied to the 
ends of the rod. We can rewrite this term in the concise form [T(s)Sq] where, 

T(s) = G x {a)2B iq + G 2 (s)2B 2 q + G 3 (s)2B 3 q. (23) 

The vector T(s) = [?i(s) T2(s) T 3 (s) T4(s)] T is the internal moment projected in 
the quaternion basis defined as a linear superposition of the internal moments 
due to elementary modes of deformation. The functionals Gi(s), G2(s) and 
G 3 (s) given by, 

G fe (s) = E k I k (/s fc (fl) - k k {s)) = E k I k (2B k q{s),q'(s) - k k {s)) (24) 

are respectively the two fiexural and torsional moments, defined as the com- 
ponents of T(s) in the local material frame. The second term in Eq. (22) is 
the work done by the operator upon a change of orientation applied along the 
rod. The elementary contribution to the integral can be rewritten J Q T(s)8qds 

where t(s) = [ti(s) T2(s) t 3 (s) T4(s)] T as a four-dimensional vector written in 
the quaternion basis that reads, 

3 

t(s) = (G k (s)4B k q'(s) + G' k {s)2B k q{s)) , (25) 
fe=i 

where G' k {s) = E k I k (2B k q' (s)q' (s) + 2B k q(s)q"(s) - k' k (s)) is the differential 
of Gfc(s) with respect to s. The quantity r(s)ds is the net moment applied on 
an infinitesimal element of the rod located between the cross-sections at s and 
s + ds. 

Before arriving to the equilibrium equations from this variation, we need to 
consider the external loads that are applied to the rod, and whose work must 
balance the variation of energy at equilibrium. Here, we consider two types of 
external loads: point forces (P(0),P(L)) and torques (M(0), M(L)) that are 
applied at the two ends s = and s = L, and distributed forces and torques 
that are applied along the length of the rod, with linear densities p(s) and m(s), 
respectively. The density of forces, p(s), can represent, for instance, the weight 
of the rod, and the density of moments, m(s), hydrostatic loadings such as the 
result of viscous stresses due to a swirling flow around the rod. The total work 
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done by these external forces upon an infinitesimal perturbation of the rod's 
configuration is, 



5W = P(0)<5r(0) + M(0)5g(0) + P{L)8r(L) + M(L)Sq(L) + 

(p(s)8r(s) +m(s)6q(s))ds, (26) 

Jo 

where Sr = [Sri & r 2 is the vector of the small arbitrary perturbations 

of the translational degrees of freedom r. According to Eq. (26), the exter- 
nal forces P(s) = [P x (s) P y (s) P z (s)] T and p(s) = [p x (s) p y (s) p z (s)] T are de- 
fined in terms of the global directions x, y, z whereas the external moments 
M(s) = [Mi(s) Af 2 (s)M 3 (s)M 4 (s)] T and m(s) = [mi(s)m 2 (s)m 3 (s)m 4 (s)] T 
are expressed in the quaternion basis i, j, k, 1 of H, defined by Eq. (8). In Sec- 
tion 5, we will show through the specific example of the writhing of a rod how 
to express physical rotational quantities (e.g. boundary conditions or external 
moments) in terms of quaternions. 

3.3. Equilibrium equations 

Thus far, we have implicitly assumed that the perturbations [5r Sq] T can 
be chosen freely. This is, however, not the case since our rod is subject to 
the kinematical constraints introduced previously in Eqs. (11) and (18). These 
constraints are imposed in the derivation of the equations of equilibrium by 
adding a number of Lagrange multipliers into the variation of the elastic energy 
£ e and external loads SW. In this Lagrangian formalism, the enforcement of 
the unicity of quaternion in Eq. (11), translates as the continuous functional 
constraint, 

C a [q(s)] = q(s)q(s) -1 = 0, (27) 

where the brackets indicate that C a depends on the function q(s), globally. 
Moreover, Eq. (18), which ensures that the directors are adapted to the Cosserat 
curve, translates to three conditions on the continuous vector-valued function, 

C„ [r(s),q(s)} = r'(s) - d 3 (q(s)) = 0. (28) 

With the expressions for the energy of an arbitrary configuration of the rod 
in Eqs. (22) and (26) in hand, the equations of equilibrium are now obtained by 
assuming that the energy is stationary under small deformations for a given set 
of boundary conditions and geometrical constraints; Eqs. (11) and (18). This is 
equivalent to requiring that the first order variation of the functionals 5£ e and 
<5W, combined linearly with the variation of the constraints SC a and SC^ over 
the interval from s = to s = L (i.e. the Lagrangian) vanish, 

S£ Lag =S£ e -SW+ a(s)SC a (s)ds + ti(s)8C n(s)ds = 0. (29) 
J a Jo 



12 



In this equation, the variation of the constraint C a {s) given in Eq. (27) takes 
the form, 

\ a(s)SC a (s)ds = / 2aqdqds, (30) 
Jo Jo 

where the scalar function a(s) is the Lagrange multiplier that imposes the norm 
of the quaternions to be one. The variation of the constraints C M ( S ) given in 
Eq. (28) reads, after integration by parts, 



/ n(s)$C-n{s)ds = [fiSr]^ - / fi'Sr + 2D(q)fiSqds 
Jo Jo 



(31) 



where the terms of the vector valued function fj,(s) = [/j, x (s) fi y (s) H z (s)] are 
the Lagrange multipliers ensuring the condition of inextensibility of the slender 
clastic rods and the operator D(q) reads, 



D(q( S )) 



93 (s) 


-94 (s) 


-Qi(s) 


94 (s) 


% 0) 


-92 (a) 


Qi(s) 


32 (s) 


-9s(s) 


92 (s) 


-91 (s) 


-94(s) 



(32) 



Now, substituting Eqs. (22), (26), (30) and (31) into the Lagrangian of Eq. (29), 
we arrive at the first variation of the geometrically constraint elastic energy of 
the slender elastic rod, 

6£ Lag = [T(s)Sq(s) + M (s)«r(a)]£ - M(0)Sq(Q) - M(L)Sq(L) 

- P(Q)8r(0) - P(L)8r(L) - [ (p(s) + fi'(s)) Sr(s)ds 

Jo 

-[ (T(s)+m{s)-2a(s)q(s) + 2D(q(s))n(s))Sq(s)ds. (33) 

The condition that the variation in Eq. (33) must vanish for an arbitrary per- 
turbations 6r(s) and Sq(s) yields the strong form of the equilibrium equations 
for our elastic rod as second-order differential equations, 

Q = p(s) + /x'(s) (34a) 
= t(s) - m{s) + 2a(s)q(s) - 2D{q(s))fi(s). (34b) 

When projected along the three directions of the global cartesian frame (x, y, z), 
the vector equation Eq. (34a) yields a set of three differential equations that can 
be interpreted as the balance of forces. The vector of Lagrange multiplier fi(s) 
measures the resultant of the contact forces transmitted through the rod's cross- 
section. Indeed, calculating the forces acting on a small element of the rod of 
length ds, we find that the element is submitted to the contact forces fi(s + ds) 
and — /x(s) from the neighboring elements, and to the external force pds. At 
equilibrium, the total forces (fi'ds + pds) is zero as described by Eq. (34a). 
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When projected along the four elements (i,j,k, 1) of the quaternion basis 
H, the vector equation Eq. (34b) yields a set of four differential equations that 
can be interpreted as the balance of moments. Working in the quaternion basis, 
it is, however, not straightforward to find an obvious physical interpretation 
for each of the terms but it suffices to say that they are related to the internal 
moments acting on a small element of the rod ds. 

For the equilibrium equations written in Eqs. (34) to be complete and well- 
posed, one must add the geometrical constraints given by Eq. (27) and Eqs. (28), 
which, in their projected and developed form, read as, 



= qf(s) + ql{s) + q£(s) + q\{s) - 1 (35a) 

= r' x (s) - 2 qi (s)q 3 ( S ) - 2q 2 ( S )q4s) (35b) 

= r' y (s) - 2q 2 ( s )q 3 (s) + 2?i («)&(«) (35c) 

= r' z (s) + qf(s) + ql(s) - q 2 3 (s) - qt(s). (35d) 



In the seven differential equilibrium equations of Eqs. (34) plus the four differ- 
ential equations in Eqs. (35), the eleven unknowns are the four Lagrange multi- 
pliers a(s), fJ, x (s), Hy(s) and fi z (s), the four rotational degrees of freedom qi(s), 
92(s), 93(s) and q^s) and the three translational degrees of freedom r x (s), r y (s) 
and t z (s). Thanks to the use of quaternions, the kinematics is geometrically- 
exact and the resultant equilibrium equations are simply polynomial since the 
highest geometric nonlinearity comes from the vector r(s) given in Eq. (23), 
which is cubic in q(s). In Section 4, while developing the numerical implemen- 
tation, we will make extensive use of this smooth and regular nonlinearity to 
efficiently compute the numerical solutions of these equations. 

So far, in Eq. (33), we have only considered the vanishing of the integral 
term. Likewise, boundary terms should also vanish since this equation is also 
to be satisfied for perturbations localized at its extremities. The first boundary 
terms, associated with rotations 6q(L) and Sq(0) yield, 

0= (T(0) + M(0))<5g(0), (36a) 
= (T(L) - M(L)) 8q(L). (36b) 

The remaining boundary terms associated with displacements dr(0) and 6r(L) 
of the ends s — and s — L, respectively, yield, 

0= QLi(0)+P(0))<5r(0), (36c) 
= {n{L) - P{L)) 8r{L). (36d) 

To provide a physical interpretation of the behavior at the boundary conditions, 
we first consider Eq. (36a). If the endpoint ,s = is free to rotate, the vector 
Sq(0) is arbitrary and one is led to the boundary condition T(0) + M(Q) = 0. 
This is the total torque applied on the section s = 0, which is the sum of the 
internal moments T(0) transmitted by the downstream part of the rod, s > 0, 
and of the moment M(0) applied by the operator. At equilibrium, the total 
torque should vanish when the end is free to rotate. If the endpoint s — is fixed, 
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the perturbations that are consistent with the kinematics are such that 6q(0) = 
and the equation is automatically satisfied. The boundary condition is then 
the one imposing the rotation of the fixed end, which leaves the total number of 
boundary conditions unchanged. The same reasoning holds for Eq. (36b) near 
the opposite end, s = L, although the total torque is now T(L) — M{L), since, 
in this case the internal moment is applied by the downstream part of the rod, 
s < L. 

The two other boundary conditions written in Eqs. (36c)-(36d) can be han- 
dled in a similar fashion. Near an end where the displacement is unconstrained, 
the total force should be zero. This total force is /n(0)+P(0) near the end s = 0, 
by a similar reasoning as above. However, the total force is (J,(L) — P(L) near 
the opposite end, s — L, given that the internal forces £t(L) are now applied by 
the downstream part of the rod, s < L. This remark validates our previous in- 
terpretation as for the physical interpretation of the Lagrange multipliers /^(s), 
fi y (s) and Hz(s); they are the internal forces along the three directions of the 
global frame that constrain the directors to be adapted to the Cosserat curve. 

Together, Eqs. (34)-(36) constitute the set of geometrically-exact cubic dif- 
ferential equations that describe the mechanical behavior of the slender elastic 
rod represented in Fig. 1(a). These nonlinear differential equations could be 
solved with classic boundary value problem algorithms upon knowing the bound- 
ary conditions in terms of external forces or kinematics. Moreover, coupled 
with traditional predictor-corrector methods, one should be able to continue, 
step-by-step, the solutions of this nonlinear elastic problem in terms of given 
geometric or mechanical control parameters [13, 17]. In the following section, in 
an alternative point of departure, we use a continuation method based on the 
Asymptotic Numerical Method (ANM) developed in the early 1990's to solve 
elastic structural problems in the early post-buckled regime [21, 22]. Taking 
advantage of the particular cubic form of the geometrically-exact equilibrium 
Eqs (34)-(36), this path-following perturbation technique will enable the deter- 
mination of semi-analytical nonlinear solution branches by inverting a simple 
stiffness matrix at each step of the continuation. This outstanding numerical 
property makes the ANM algorithm highly robust and computationally efficient 
at determining the various equilibria of our slender elastic rod. 

4. Numerical method 

In this Section, we solve the differential equilibrium equations Eqs. (34)- (36) 
using a finite element-based semi-analytical path- following method. We first 
approximate the continuous degrees of freedom using finite differences approx- 
imation to interpolate the mechanical and geometrical variables at each nodes 
and elements. Thanks to the quaternion formalism introduced above, the equi- 
librium equations can be reduced to an algebraic set of quadratic equations 
by considering the flexural and torsional internal moments as unknowns. This 
quadratic form is particularly well suited to ANM which is a semi-analytical 
continuation algorithm to compute the branches of solution of a set of nonlinear 
polynomial equations. To follow all the bifurcated branches, we show how the 
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local stability of the computed equilibria can be assessed by using the second 
order conditions of constrained minimization problems. Finally, we describe the 
implementation of our algorithm into MANlab, a free and interactive bifurcation 
analysis software based in MATLAB. 



4--1- Discretization 

In order to compute the equilibrium equations, Eqs (34)- (36), we first explain 
how to discretize the main function unknowns such as strain rate vector k(s), 
material frame (di(s) d 2 (s) ds(s)), positional and rotational degrees of freedom 
r(s) and q(s) or Lagrange multipliers a(s) and n(s). 

The position of the rod is represented by discretizing its centerline into N 
elements separated by N + 1 spatial control points, r(s l ) = r % = \r % x r l y r\\ in 
1Z 3 , located by the discrete curvilinear coordinate s l as illustrated in Fig. 3(a). 
The spatial derivative of the positional degrees of freedom is approximated by 
the forward finite difference between two successive nodes, 

Q i\ r i+l _ r i 

rV) « ' "~ ' Z' ' V } = (37) 



(s l + ds l ) - r(s l ) _ 

ds l ds i 

where ds l = ||r ,+1 — r' l \\ = L/N is the length of the i th element and L is the 
total length of the rod. To ensure the inextensibility condition of our rods, ds l 
is constant upon deformation and the stretch along the centerline is forced to 
verify \\r' {s l )\\ = 1. 

The orientations of the centerline elements arc represented by employing N 



Ti <?2 % ll 



is the 



material frames R(q-') = d\dl l d\ in lZ 3x3 where g J = 
set of quaternions associated with each element j. According to Eq. (12), the 
directors of the j th element are vectors in 1Z 3 represented at the midpoints on 
the centerline segments [see Fig. 3(a)] such that, 



d{ 



1 1 



1 7 

Q2I2 



4,4, 



1 1 



^(qW 2 + ikl) 

L 2(qiq> 3 - q J 2 qi) 



2(44 - 44) 
-44 + 1W2 - 
2(44+4,4) 



1 7 

q\q\ 



^ '/)'/:; + ikl) 

2(<f 2 4 - 44) 



1 1 

-qWi 



1 1 



1 7 



7 7 

qiqi 



(38) 



Replacing the quaternions functions q(s) by their discrete counterparts q 3 in the 
expression of strain rates given in Eq. (16), we can write the discrete material 
curvatures n\ , k j 2 and the twist k j 3 expressing the extent of rotation around the 
directors, d\, d 2 and d 3 , between two successive elements [see Fig. 3(b)] in the 
form 

K> k = 2B k q 3 q' 3 for k=l, 2, 3. (39) 

In Eq. (39), we introduced the average and the spatial derivative of the rotational 
degrees of freedom of the j th element q 3 as, 



7J+ 1 



q> = 



(40) 



16 




bi) b2) b3) 




Figure 3: Finite element discretization method, (a) The centerline of the rod is discretized 
into N elements separated by N + 1 nodes, r J . We also consider N material frames H(q J ) 
to express the orientation of the j th element, (b) The change of orientation between two 
successive elements j and j + 1 is expressed by the N — 1 discrete material curvatures at the 
interconnected nodes: (bl) rej around the director (dj +d| +1 )/2, (b2) k\ around the director 
(dj + df 2 +1 )/2, (b3) k{ around the director {d{ + d{ +1 )/2. 



„i+2 



„i+2 



= L/N, taking into consideration 



where ds 3 = l/2(||r' _r - i -r' Ti || + ||r' T - i -r 
the rod's inextensibility condition. 

In a similar fashion, replacing the continuous function q(s) and its deriva- 
tive q'(s) by their discretized counterparts, q 3 and q' 3 , respectively, given in 
Eq. (40), the first variation of elastic energy previously given in Eq. (21) can be 
approximated by a Riemann sum over the elements from j = 1 to N, 



3 N-l 



2. 2. 2G i 



(B k q j 6q j+1 - B k q 3+1 8q 3 ) ds j . 



(41) 



k=l j = l 



In this equation, the N vectors Sq 3 — Sq{ Sq? Sq^ Sq 3 4 are the discrete version 

of the perturbed rotational degrees of freedom Sq(s) and are associated with 
each element j. The 3(JV — 1) constants Gi are an approximation of the flexural 
and torsional internal torques Gfc(s) given in Eq. (24), which read 



G{ = E k I k [ B k (q> 



ds 3 



(42) 



and arc defined between two successive elements for 1 < j < N. 
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Replacing q(s) and r(s) by their discrete counterparts at each element qi and 
each node r l , respectively, the variation of the work done by the external forces 
and torques previously given in Eq. (26), can be approximate by its discrete 



SW « P^dr 1 + M°Sq 1 + P L 6r N+1 + M L Sq N + 

N JV-1 

^p'<5rW + Y rn 1 8q i ds j , (43) 

i=2 j=2 

where Sr l = [Sr l x 5r y 8r\\ T is the vector of the perturbed displacement at each 
node. In this equation, we have introduced the vector of point forces at the two 
ends P° = [P° P° P°] T and P L = [p£ P£ Pf] T and the vector of density of 

external forces at each node p 1 = [p* p y p],] for 2 < i < N defined in term of 

the global directions x, y, z. In the discrete version of 8W, Eq. (43), we have 

also introduced the torques applied at the two ends M° = [M? M§ M| Mf\ T 

and M L — [Mjp M% M±\ T and the vector of density of external moment 
r --\T 

at each element m? = mj TO3 to^ for 2 < j < N — 1, also expressed in 
the quaternion basis. 

Before we derive the algebraic system of equilibrium equations, we still need 
to write the discrete form of the variation of work due to the geometrical con- 
straints in Eqs. (27) and (28). Replacing q(s) by its discrete counterpart, qi , 
we can expand Eq. (30) in the form of a Riemann sum, 

rL n 

/ a{s)SC a (s)ds « V2aV<5qW, (44) 
Jo i=1 

where a? is a discrete scalar at each element j, approximating the continuous 
Lagrange parameter given in Eq. (30). Introducing the N vectors of Lagrange 

parameters /i J = \jj? x H J y nl\ T which prescribe that each element is parallel to 
the tangent r'(s l ) given in (37), we can rewrite Eq. (31) in its discrete form, 



Jo 



N N 

Y (m* - A 4 * 1 ) Sr l + Y O q .frSq d.s 1 . (45) 

i=2 j=l 

where the vector of Lagrange parameters, or internal contact forces, jtx(s), has 
been approximated by the backward finite difference between each successive 
elements, 
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and D(q : ') is the discrete counterpart of D(q(s)) introduced in Eq. (32) which 
reads, at each element, 
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(47) 



As we did previously in the continuum case, we now require that the discrete 
variation of the Lagrangian 5£Lag (given in Eq. (29) as the sum of Eqs. (41), 
(43), (44) and (45)) vanishes for any arbitrary perturbations 6r l and 8qi . This 
condition yields the set of algebraic equilibrium equations of the discrete un- 
shearable and inextensible slender elastic rod. The condition that this variation 
is zero for any perturbed displacements Sr l leads to the balance of forces as a 
set of algebraic equations, 



p l L/N + ^ 



.2-1 



for i = 1 
/x " ^ for 2 < i < N 
(j~ 1 =0 fori = JV+l. 



(48a) 
(48b) 
(48c) 



Projected along the three directions of the global cartesian frame (x, y, z), 
Eqs. (48) yield 3(7V + 1) linear equations for the 3iV unknowns \i d . In the 
limit N very large, these equations converge to the continuous differential equa- 
tions Eq. (34a). The condition that the variation of the Lagrangian is zero for 
any arbitrary perturbations 8q J leads to the balance of moments as a set of 
discrete algebraic equations, 



tJ 



■Af 



M 



2D{q j )^ 
2D(q>)^ 
2D(q j )fj, j 



2a V = for j = 1 (49a) 
2a V =0 for 2 < j < N - 1 (49b) 
2aV = for j = N. (49c) 



In these equations, t j 



•i '2 '3 '4 



is the vector of net internal moment 



applied on the element j written in the quaternion basis, 

rj = ELi G i 2B ^ J+1 for i = 1 ( 49d ) 

t° = {G{7 1 2B k q^ 1 - G j k 2B k q ]+1 ^j for 2 < j < N - 1 (49e) 

t 3 = _ G{7 1 2B k q 3 - 1 for j = N. (49f) 

and 2D(q 3 )/j, : ' is the moment resultant from the internal contact forces ^ 
applied on the element j. In the limit of large N, Eqs. (49) converge to the 
continuous differential equation Eq. (34a). Projected along the four elements 
k, 1) of the quaternion basis H, Eqs. (49) yield 47V nonlinear equations for 
the 47V unknowns q J and N unknowns a J . The missing equations required to 
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compute all the unknowns are given by the geometrical constraints Eqs. (35) 
which can be rewritten in the algebraic form, 



Wif + + iti? + (d) 2 - 1 - for 1 < j < N 
r i+i _ r 3 _ dids j = for 1 < j < N, 



(50a) 
(50b) 



where we used the forward finite difference in Eq. (37) to approximate r'(s). 

The 47V geometrical constraints given in Eq. (50), together with the 77V + 3 
equilibrium equations Eqs. (48) and (49) form the set of algebraic equations 
describing the constrained equilibrium configuration of the rod represented by 
the 7N + 3 degrees of freedom and and the AN Lagrange parameters a- 7 
and /x- 7 . 

We highlight the fact that the only approximations made in the above equa- 
tions arise from the finite element discretization since the initial continuous for- 
mulation is geometrically-exact due to the use of quaternions. Furthermore, it is 
remarkable to notice that the equilibrium configurations of the extremely twisted 
and bended elastic rod can be represented by the smooth polynomial equations 
Eqs. (48)-(50). In the next section, we exploit the particularly smooth nonlin- 
earities of the equilibrium equations by using Asymptotic Numerical Methods 
(ANM) [21, 22, 23] which are efficient path-following techniques that give access 
to semi-analytical solution branches of polynomial nonlinear algebraic systems. 

4-2. Asymptotic Numerical Method 

We now explain and adapt the particular ANM introduced in [22] for solving 
the equilibrium equations of slender elastic rods described above. This ANM 
is a perturbation technique allowing for the computation of a large part of a 
solution branch of quadratic algebraic system of equations with only one stiff- 
ness inversion. Applied in a step-by-step manner, one can compute a complex 
nonlinear branch by a succession of local asymptotic expansions and thus de- 
termine a semi-analytical bifurcation diagram. Because of the local analytical 
representation of the branch within each step, this continuation technique has a 
number of important advantages when compared to classical predictor-corrector 
schemes [22]. In particular, the algorithm is fully automatic, remarkably robust, 
and faster than incremental-iterative methods. 

To apply the asymptotic numerical method to the mechanics of elastic rods, 
we first rewrite the algebraic nonlinear systems of equilibrium equations Eqs. (48)- 
(50) in the compact form, 



number of elements of the discretized rod, A is a scalar control parameter (usu- 
ally a mechanical or geometrical parameter of the physical problem such as the 
rotation or displacement at one end of the rod) and w is the vector of unknowns 
which, in our case reads, 




w = [r 1 q 1 . . . 




11 N N1 T 

a ft ... a fj, \ . 



(52) 
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According to the discretization presented in Section 4.1, w is a vector of size 
UN + 3 which includes the 7N + 3 mechanical degrees of freedom separated 
into positions r % and quaternions g 3 '. Moreover, the 47V Lagrange parameters 
are required to impose the geometrical constraints. 

In what follows, in a process that we refer as recasting, we now transform 
Eq. (51) into a quadratic form, which is a particular framework of the ANM that 
allows us to formally and systematically write a large class of physical problems 
including rods [22, 23]. Given the original cubic form of Eq. (51), this quadratic 
recast is achieved introducing a new vector of unknowns u of size 14 A + 1, 

u=[wG 1 1 G 1 2 G 1 3 ...G 1 N G 2 N G 3 N \] T , (53) 

which includes the initial vector of unknowns w given in Eq. (52), the control 
parameter A and where we added the 3(A — 1) fiexural and torsional internal 
torques G°Aqi) introduced in Eq. (42). Using the new vector u instead of w, we 
can recast the cubic nonlinear vector valued function / (w, A) given in Eq. (51) 
into the quadratic form, 

/ («) =L + L{u) + Q(u, u) = 0, (54) 

where / (it) is a vector in 1Z 14N since we added the 3(A — 1) nonlinear quadratic 
equations Eq. (42) to f (w), Lq is a constant vector and L(m) and Q(»,») are 
a linear and bilinear vector valued operators, respectively. The expression of 
/ (if) , representing the equilibrium equations of our inextensible and unshear- 
able elastic rod given in Eqs. (48)-(50), is provided in Appendix A. In Section 5 
where we will apply our method to the quasi-static writhing of a double-clamped 
elastic rod, we will illustrate, by way of example, how to introduce the boundary 
conditions and the control parameter in the vector / (u) of Eq. (54). 

We can now proceed and compute the solutions u of the set of quadratic 
equations Eqs. (54) with the asymptotic numerical method. This technique is 
based on the perturbation of the vector of unknowns u in terms of a path- 
parameter a in the form of the asymptotic expansion, 

m 

u(a) — u + a p u p , (55) 
p=i 

where Uq is the starting fixed point, solution of Eq. (54), m is the truncation 
order of the power series and a is the path-parameter which will be formally 
defined below. Replacing u(a) by its asymptotic expansion Eq. (55) in the 
quadratic form Eq. (54), we obtain the quadratic Taylor series in the neighbor- 
hood of Uq, 

f (u{a)) = L + L(u ) + Q(u 0l u ) 

+ a[L(u l ) + 2Q(u ,u l )] (56) 

Em •r-^p—l 
a p L(u p ) + 2Q(ii , Up) + > Q(iti,ii p _ 2 ) 
p— 1 * J i—1 
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Recalling that it is a solution of Eq. (54), we can rewrite Eq. (56) in the 
form of a power series of a quadratic vector valued function f p of size 147V, 

/ («(o)) = o/i + a 2 / 2 + . . . + a m / m = 0. (57) 

Since Eq. (57) has to be verified for every value of a, we need f p = for every 
order p < m. This leads to to linear systems in u p in the form, 



Vp e [l...m], 



where, due to the particular quadratic form of Eq. (54), the Jacobian matrix of 
f(u) evaluated at the initial solution vector u Q reads, 



df 
du 



Up = L(u p ) +2Q(u 0) u p ) ) (58b) 



and the nonlinear vector f™ 1 on the right-hand side of Eq. (58a) consists of a 
quadratic sum that only depend on the previous order, 

/£' =0 for p = 1 (58c) 

fp 1 = -Y? 1 Q{Ui,u p -i) for 1< p < to. (58d) 

^ * — ^ i— 1 

The original nonlinear problem in Eq. (54) has thereby been reduced to a definite 
set of m linear systems given in Eq. (58a) where the matrix on the left-hand 
side is identical for each order. 

However, each linear system in Eqs. (58) is, so far, under-determined since 
the dimension of / is 14iV whereas the size of the vector of state variable at 
order p, u p , is 14 N + 1. The remaining equation is provided by the definition 
of the path parameter a as defined in [22] . We consider a measure that includes 
the entire set of physical unknowns and that is also robust towards limit and 
bifurcation points, i.e. an arc-length measure. Mathematically, we identify the 
path parameter a as the projection of the vector of state variables increment 
u — Uq on the normalized tangent vector U\ [see Fig. 4(a)], 

a = (u(a) — u ) T Ui. (59) 

Replacing u(a) by its asymptotic expansion, Eq. (55), in Eq. (59), we obtain, 

au\ T Ui — a + a 2 Ui T U2 + . . . + a m Ui T u m = 0. (60) 

Verifying Eq. (60) at every power of a provides us with the supplementary 
equations at every p, 

Ui T u p — S p i for 1 < p < to, (61) 

where S p ± is the Kronecker delta, 5u = 1 at the first order p — 1 and it is zero 
otherwise. 
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Figure 4: (a) The path parameter a is identified as the projection of the solution branch on the 
tangent vector. This projection is normalized by the length of the tangent vector [u(l) A(l)] 
which is set to unity, (b) Behavior of a power series close to the radius of convergence. The 
considered nonlinear unidimensional equation is f (u(a)) = u(l + a) — 1 = 0. Its unique 
solution u(a) = l/(l + a), represented as solid line, can be expanded asymptotically as u(a) RJ 
1 - a + a 2 + . . . + (-l) m a m + (-l) m+1 a m+1 . The asymptotic expansions for m = 5, 10, 15 
and 20 are represented as different lines. We see that u(a) can be approximated up to the 
radius of convergence a = 1. Applying our ANM method, the step length calculated through 
Eqs. (64)-(66) would give a max = eVCm+l) = Q 4642 with e = lx w -7 and m = 2 Q. 



Finally, the original nonlinear problem in Eq. (54) has now been transformed 
in the well-posed m linear systems in 1Z 14N+1 , 
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for 1 < p < m, 



(62) 



with a unique solution u p that we can solve iteratively since each of these vectors 
is defined with the solution of the previous order according to the definition of 
given in Eq. (58). In this linearized numerical problem, the only matrix to 
inverse is the one on the left-hand side of Eq. (62) since it is the same at every 
order p. This is in striking contrast with classical predictor-corrector methods 
where one needs to actualize the Jacobian for every linear systems [13]. 

For practical purposes, one will inverse the matrix and compute the unknown 
Ui at first order independently and then compute the higher orders u p for 
1 < p < m from the well-defined systems in Eq. (62) . 

Once each u p has been found, we still have to estimate the validity domain 
of the asymptotic expansion since Eq. (54) can only be true for values of the 
perturbation parameters a inside the radius of convergence of the power series 
given in Eq. (55) [see Fig. 4(b)]. A simple, robust and accurate way of calculating 
an approximation of the convergence radius a max , explained in detail in [22], 
is to assume that a solution branch is acceptable as long as the norm of the 
nonlinear (14iV + l)-dimensional vector field / (u(a)) is less than a tolerance 
criterion e, 

Vae[0a max ], ||/(o)||<e, (63) 
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where e determines the accuracy of our numerical results. We have computed 
the Up according to the power series expansion of / (a) given in Eq. (57) so that 
the norm of / is zero up to the truncation order m. Consequently, the residue 
of this series is given by the norm of / (a) for p > m. Assuming that the order 
m + 1 dominates in the residue, we obtain the relation between the norm of 
/ (a) and the vector at the order m + 1, 

||/ (a) || * a m+1 \\f m+1 \\ (64) 

where f m+1 = f£| Mo «m+i-/m+i according to Eqs.(57)-(58). Replacing ||/ (a) || 
by its definition Eq. (64) in Eq. (63) leads to, 



||/(a)||<e a <[yT^) ■ ( 65 ) 

which sets an upper limit to the path-parameter a. While truncating the asymp- 
totic series u{a) given in Eq. (55) at the order m, we implicitly assumed that 
u m +i = so that f m+ i = Sm+\- Replacing f m+1 by the nonlinear term f%f+i 
in Eq. (65), we obtain an estimation of the maximum step length, 

i 

(66) 



m+1 1 

For practical purposes, when applying Eq. (66) to any quadratic vector valued 
function f(u(aj), we found that / (u(a max )) rj e. In general, the power series 
in Eq. (55) converges slowly, close to the radius of convergence [23] . Decreasing 
e leads to a diminishing of a max but, more importantly, results in an increased 
accuracy of the computed series. An optimal value of m and e in terms of 
convergence of the asymptotic series, accuracy of the solution given by Eq. (63) 
and size of a max is found empirically for the following range of parameters [23] : 
e = 1 x KT 7 and 15 < m < 20. 

The power series expansion given in Eq. (55) and computed with the m linear 
systems in Eq. (62), together with the maximum step size a max given by Eq. (66) 
define a portion of the nonlinear equilibrium branches of the slender elastic rod 
in terms of a given control parameter A. The next step of our calculation, the 
continuation of the solution branch, is now computed by applying the present 
asymptotic numerical method taking u(a max ) as the new starting equilibrium 
u of the new portion. A complete solution branch is therefore constructed as 
a succession of semi-analytical portions in the form of Eq. (55), whose length is 
automatically determined through the estimation of the convergence radius of 
each power series as sketched in Fig. 5(a). Unlike classical predictor-corrector 
methods [15, 17], our step length is adaptive; it is naturally large for weakly 
nonlinear solutions and becomes shorter when strong nonlinearities occur. As 
a consequence, the automatization of the continuation method is significantly 
easier and more robust than with standard predictor-corrector methods. 

One would expect that the residue ||/ (u(a)) |j would increase progressively at 
every continuation step so that the accuracy of the new starting equilibrium Uq 
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Figure 5: (a) Schematic five continuation steps with the asymptotic numerical method. The 
step length of each portion of branch is determined a posteriori by analyzing the validity of 
the asymptotic solution. The highly nonlinear solution is known analytically for each portion, 
(b) Branch switching through the perturbation method given in Eq. (67). Changing the sign 
of the intensity of perturbation, c, allows us to explore all the branches after a bifurcation 
point. 

would gradually decrease [23]. In practice however, it is rare to see the residue 
increase up to lOe, especially given the smooth nature of the nonlinearities of 
our equilibrium equations Eqs. (48)-(50) for thin elastic rods. In Section 5, 
where we implement this continuation method to a series of specific test-case 
problems, all the bifurcation diagrams are computed with a residue smaller 
than e = 1 x 10~ 7 , with no correction step (note that a correction step may 
be necessary in the general ANM framework in the case of non-polynomial 
nonlinearities [44]). 

When the accuracy e is set to be small enough by the user, the ANM is 
able to follow the branch whenever bifurcation are encountered [45]. This is 
a remarkably robust property for a path-following algorithm, especially when 
compared to predictor-corrector techniques which typically would systematically 
bifurcate because of the discrete nature of their continuation steps. Nevertheless, 
we now need a special procedure to switch branches in order to determine the 
full bifurcation diagram. A classic strategy is to slightly modify the original 
equilibrium equations Eq. (54) by adding a low-norm perturbation vector, 



where fp (it) is the perturbed problem, P a normalized vector of constant 
random numbers and c is the intensity of the perturbation. This additional 
perturbation procedure transforms the exact bifurcation into a perturbed bi- 
furcation [see Fig. 5(b)]. The idea is to use the perturbed branch to bifurcate 
on the non-crossing branch [46]. Changing the sign of c allows us to explore 
a symmetrical quasi-bifurcation as represented in Fig. 5(b) with the perturbed 
branches P = +c and P = —c. Finally, in order to transition from the original 



f P (u) =L + L(u) + Q(u, u) + cP, 



(67) 
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to the perturbed problem, or vice versa, a correction step is mandatory (any 
predictor-corrector methods would be efficient since the perturbed solutions are 
very close to exact ones). The combination of different (positive and negative) 
values of the intensity of the perturbation and several correction steps allows 
one to explore the full bifurcation diagram of the slender elastic rod described 
by the equilibrium equations Eqs. (48)-(50). 

So far, we have presented the ANM method in the context of thin elastic 
rods. We can now compute the bifurcation diagrams of our slender elastic rod 
under various mechanical and geometrical environments. However, the final 
crucial step of determining the stability of the solution branches is still missing, 
which is the focus of the following section. 

4-. 3. Stability analysis 

Determining the local stability of equilibrium branches is crucial for the 
physical understanding of the mechanical behavior of slender elastic rods, one 
of the main motivations being that locally unstable branches cannot be observed 
experimentally, and must therefore be classified. Another advantage for gaining 
knowledge on the stability of a solution is that the loss of local stability is 
often associated with a bifurcation point. Assessing the stability is then useful 
to detect bifurcation points and navigate through the bifurcation diagram as 
illustrated in Fig. 5(b). 

The equilibrium equations on their own are not sufficient to determine the 
local stability of the solutions; we also need to compute wether the solution is 
a local minimum or maximum of the system's energy. For practical purposes, 
we need to derive the second-order conditions of the geometrically constrained 
energy; a theoretical and numerical procedure that is well established [14, 47]. 

In the previous sections, we have shown how to compute and follow the 
branches of solutions of the nonlinear algebraic equilibrium equations given in 
Eqs. (48)-(50). Since the resulting computed bifurcation diagrams are semi- 
analytical, we are therefore able to evaluate the solution u e given in Eq. (53) at 
any finite value of the control parameter A e . Let x e , a e and fi e be the vectors 
of degrees of freedom and Lagrange multipliers respectively, associated with the 
vector solution u e , such that, 



The vector X(3 IS £1. solution of the equilibrium equations Eqs. (48)-(49) for A = A e 
that satisfies the functional geometrical constraints given in Eq. (50) . Rewriting 
Eqs. (48)- (50) in an energy minimization framework, x e is the actual solution 
of the n — 7N + 3-dimensional constraint minimization problem, 




U\...r N x r» rf q? q» g» qf r£™ r^f , (68a) 



(68b) 



N 3N 



V (£ e (x e ) + W (x e )) + £ oc\VC a + £ /4VC£ = 0, 



(69a) 



i=i j=i 



26 



subject to the m = 47V functional constraints, 
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for the positional nodes 1 < j < N. In Eq. (69), V is the gradient operator 1 , 
the real- valued function £ e (x e ) measures the amount of elastic energy stored in 
the rod at equilibrium, 

s e (x e ) = J2 £ ( B » («* + I W +1 *£) 2 dsJ . (™) 

fe=l j=l ^ ' 

and W (x e ) quantifies the total work of external forces and moments, 
W (x e ) = PV + MV + P^^ 1 + M^q^-I- 

N N-l 

y^ j p i r i ds + ^ rn 3 q j ds. (71) 

Finally, the vectors C Q = [C* C 2 a . . . C£] T and C M = [C^ C^... C% C% C^ z ] T 
in Eq. (69a) represent the geometrical constraints that ensure the norm of 
quaternions to be one and the inextensibility of the oriented Cosserat rod re- 
spectively, and should be very close to zero at equilibrium. 

According to the necessary and sufficient first order conditions of constrained 
minimization problems [14, 47], the vector solution x e is a local extremum (a 
minimum or maximum) of the total energy £(x e ) — £ e (x e ) + W(x e ) subject to 
the m constraints in Eqs. (69b)- (69c). Supposing also that the n x n matrix, 

N 3N 

L(x e ) = V 2 (8 e (x e ) + W (x e )) + J2 ^V 2 C' a + ]T y? : V 2 C£, (72) 

i=i j=\ 

where V 2 is the Hessian operator 2 , is positive definite on the m-dimensional 
subspace M = {y : X?h(x e )y = 0} with h(x e ) — [CaC^, that is, for y e M 



x For a real-valued function / £ C' 1 on TV 1 such that f(x) = f(x\ ), we define 

the gradient of / to be the n-dimensional vector, 

T 

It I I . I" I ( / j I . » ' J It' !'.*"! 



df(x) df{x) df(x) 
dxi 8x2 dx„ 



2 We define the Hessian of / at X (/ 6 C 2 ) to be the n X n symmetric matrix denoted 
V 2 /^), 

'd 2 f{xY 



V 2 /» 



dxidxi 
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and t/ / that holds y T L(x e )y > 0, then, according to the second-order 
necessary and sufficient conditions, strict local minimum of £ (x e ) subject 

to C a = and C M = 0. 

The matrix L{x e ) is the matrix of second partial derivatives, with respect 
to x, of the discrete counterpart of the Lagrangian given in Eq. (29). When 
restricted to the subspace M that is tangent to the constraint surface and which 
we denote by Lm, L{x e ) plays the role in second-order conditions directly 
analogous to that of the Hessian of the objective function in the unconstrained 
case [14, 47]. The eigenvalues, a i: and associated eigenvectors, yt, of Lm, 
determine the local stability of the solutions of the constrained minimization 
problem. Mathematically, Lm is a (n — to) x [n — to) matrix defined, at each 
equilibrium point 

L M (x e ) = (ker(Vh)) T L (ker(Vfc)) , (73) 

where ker(») denotes the kernel operator. Analyzing the n — m eigenvalues <Ji 
gives us information on the behavior of the associated perturbation yi (t) in the 
neighborhood of the equilibrium x e . According to Lyapunov's theorem [48, 49]: 

• If <Ji > for all £ € [1 . . . n — m], all the perturbations vanish, y(t) — > 0, 
when t — > +oo, and the equilibrium is locally asymptotically stable. 

• If one index i G [1 . . . n — m] exists, for which Oi < 0, one perturbation 
diverges, y(t) — > +oo when t — > +oo, and the equilibrium is locally unsta- 
ble. 

• If <Ji > for all i £ [1 ... n — to] and if there exists one index k such 
that Ofc = 0, the first order is insufficient to draw conclusions on the local 
stability of the equilibrium. In that (.'MSG , M perturbation at higher order 
is necessary. 

Applying this method to a sufficient number of fixed points x e along the 
equilibrium branches, computed with the previous ANM method, allows us to 
determine the stability of the bifurcation diagram. The previous finite ele- 
ment discretization presented in Section 4.1, together with the ANM algorithm 
described in Section 4.2 and the previous stability method, completes the semi- 
analytical continuation technique that we developed to compute and follow the 
equilibrium branches and the stability of an inextensible slender elastic rod 
undergoing extreme displacements and rotations. The combination of the con- 
ciseness and relative simplicity of our method offer the opportunity for it to be 
implemented in any programming language. In the following, we briefly present 
MANlab [26, 27], an open-source bifurcation analysis software that provides a 
convenient framework to implement the previous numerical methods. 

4-. 4- MANlab: an open-source bifurcation analysis software 

MANlab is an interactive software package for the continuation and bifurca- 
tion analysis of algebraic systems, based on ANM continuation, and first released 
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in 2009 [26]. Thanks to the implementation of most of the ANM equations in 
MATLAB using an object-oriented approach [27], MANlab makes it simpler for 
the user to solve the system of Eqs. (48)-(50) and the stability of the solutions 
given by the second-order condition through Eq. (72). MANlab has a graphical 
user interface (GUI) with buttons, on-line inputs and graphical windows for 
generating, displaying and analyzing the bifurcation diagram and the solutions 
of the system. A unique identifying feature, when compared with other contin- 
uation codes, is that its computational efficiency, highlighted above, allows for 
interactive control of the continuation process. The full interactive and semi- 
automatic procedure consists of computation of a portion of a branch, choice 
of a new branch at a bifurcation point, reverse direction of continuation on the 
same branch, jump capability between solutions, visualization of user-defined 
quantities at a particular solution point, selection and deletion of a branch, 
or of one of its portion, possibility of correction step with a Newton-Raphson 
method and determination of the local stability of the solution. 

To enter the system of equations, the user simply has to provide the three 
vector valued Matlab functions corresponding to the constant, linear and quadratic 
operators Lq, L(u), and Q(u,u) given in Eq. (54). To assess the local stability 
at each computed solution point in MANlab [50] , one can also provide the Hes- 
sian of the constrained Lagrangian restricted to M, Lm(x e ) given in Eq. (73) 
and the package will automatically compute the eigenvalues of the linearized 
problem according to the previous section. Thanks to the flexibility offered 
by the MATLAB environment, users become rapidly familiar with MANlab. 
Calling of external routines such as finite elements codes is also possible. 

In the following section, we validate our semi-analytical continuation method 
by using the MANlab package to simulate a precision model experiment; the 
quasi-static writhing of a double-clamped slender elastic rod, which equilibria 
are solutions of the discrete equilibrium equations given in Eqs. (48)-(50). 

5. Following the equilibria of an extremely twisted elastic rod 

Having introduced the general theoretical and numerical framework to con- 
tinue the equilibria and stability of slender elastic rods, we proceed by imple- 
menting the specific problem of the writhing (extreme twisting) of a clamped 
clastic rod. Even though this fundamental problem appears seemingly simple, it 
can display an array of complex behavior with intricate bifurcation diagrams and 
has received significant attention in the literature [18, 51, 29, 30, 52, 53, 28]. 
The writhing of an elastic rod is therefore an ideal scenario to challenge our 
theoretical and computational framework by contrasting the numerical results 
with our own precision model experiments that were especially developed for 
the testing and validation of our continuation method. 

In this Section, we first present our apparatus and model experiments which 
consist of quasi-statically increasing the rotation angle at one end of a slender 
elastic rod fixed between two concentrically aligned horizontal clamps. One of 
the originalities of our experiments is that we fabricate our own elastic rods, 
enabling us to accurately target their material and geometrical properties. In 
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particular, we have full control in setting their intrinsic natural curvatures. After 
describing how to account for the kinematic boundary conditions and control 
parameter specific to this writhing problem in the numerical model described 
in Section 4, we compare a series of experimental and numerical results for two 
different elastic rods: a straight rod with no natural curvature and a curved 
rod. 



5.1. Manufacturing of rods and experimental apparatus 

Our rods are cast by injecting vinylpolysiloxane (VPS), a two-part silicone- 
based elastomer, into a flexible PVC tube of inner and outer diameters D] = 3.1 
mm and Do — 5 mm, respectively. The PVC mold is first wound around 
a cylinder of external radius R e and then injected with VPS, which eventually 
cross-links at room temperature [see inset of Fig. 6] . After a setting period of 24 
hours, to ensure complete curing of the polymer, the outer flexible PVC pipe is 
cut to release the inner slender VPS elastic rod with a constant natural curvature 
ki(s) = l/(R e + Di/2) and a circular cross-section R = Dj/2 = 1.55 mm. The 
rod's second moments of area are I\ = I-x = irR 4 /4 and 1% = J = ttR 4 /2. We 
measure the Young's modulus of the elastomer to be E = 1000 KPa, a volumic 
mass p = 1200 kg/m 3 and a Poisson ratio of v « 0.5, so that its shear modulus 
is G = E/2(l + v) = 305 KPa. 





Figure 6: The writhing experiment. A L = 30 cm long elastic rod with a circular cross section 
of radius 1.55 mm, Young's modulus E = 1000 KPa, volumic mass p = 1200 kg/m 3 and a 
custom made constant natural curvature k\ is fixed at both ends between two concentrically 
aligned drill chucks separated by a distance d = 22 cm. Our experiment consists of quasi- 
statically increasing the rotation $ at one end and investigating the evolution of equilibrium 
state with the control parameter Inset courtesy of Khalid Jawed: Fabrication process of an 
elastomeric rod. During working time, the PVC tubes containing the silicone-based elastomer 
are wound around cylinders of external radius R e and left in that position for 24 hours. After 
demolding, it confers to the rod a constant natural curvature k\(s) l/R e - 
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The cast rod (L — 30 cm long) is then attached between two horizontal 
concentric drill chucks of a lathe, separated by a distance d — 22 cm. A photo- 
graph of the side view of the experiment is presented in Fig. 6. The boundary 
conditions of the rod are set to be rigidly clamped at both ends. For future 
representation of the rod configurations, we choose the origin of the cartesian 
frame (x, y, z) to be located at the clamp at the left extremity of the rod [see 
Fig. 6]. The clamp located at the origin, at the curvilinear coordinate ,s = 0, 
is completely fixed but the other clamp, located at s = L, can be rotated with 
respect to the y-axis, thereby imposing a rotation angle $ [see Fig. 6]. Initially, 
for $ = 0°, we ensure that the sign of the intrinsic curvature ki is such, that 
the rod naturally bends downwards, in the direction of gravity and that the 
difference between twist angles, ^3(5), at both ends of the rod, is zero. In that 
configuration, the equilibrium shape of the clamped rod is close to a planar 
inflectional clastica as theoretically described in [54]; the only difference arising 
from the effects due to gravity which induces a catenary-like configuration. Our 
writhing experimental protocol then consists of quasi-statically increasing the 
rotation angle, $, at s = L and quantifying the evolution of equilibrium states 
as a function of this control parameter, $. A variety of measurements on the 
configurations of the rod are performed by imaging the top of the experiment 
(using a Nikon D90 SLR camera) and subsequent image processing. 

5.2. Modeling of the boundary conditions for the writhing configuration 

Before we can proceed with a direct comparison between experimental and 
numerical results, we first need to precise how to account for the specific kine- 
matic boundary conditions and control parameter, relevant to this specific 
writhing configuration, in our general numerical framework presented in Section 
4. This specific implementation will serve as an example, which, following the 
series of procedures and rationale described below, can be extended to other 
kinematic conditions to solve a variety of other problems involving thin rods. 

Representing the slender elastic rod of Fig. 6 by the discrete 3D Cosserat 
curve of Fig. 3 and applying the numerical method of Section 4, we can write 
its equilibrium equations in the quadratic form / (it) defined in Appendix A. 
In the writhing experiment, gravity is the only external force applied to the 
rod. This gravitational force is represented by the weight of each element, 
reported at each node. In f (u) given in Eq. (A.l), we can therefore write 
p o = pL = -i/2pg 7 rR 2 L 2 /N 2 e z at both end nodes and p l = -pgnR 2 L/Ne z 
for all the internal nodes 1 < i < N+l, where g = 9.81 m/s 2 is the gravitational 
acceleration, N is the number of segments and e z is the unit vector in the z- 
direction. Since there are no external moments, we can also write m 3 = for 
all the internal elements 1 < j < N and M° = M L = at both ends. 

In addition to the mechanical parameters which are general to thin rods, 
we also need to account for the kinematic boundary conditions and control 
parameter $, specific to the writhing experiment, which are not included in the 
general formulation of / (it) given in Eq. (A.l). At the left extremity of the rod 
(s = 0), the first node must be fixed to the origin and the first element, which 
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Figure 7: Schematics of the rotational boundary conditions of our writhing experiment, (a) 
At s = 0, we need to orient the first element along the {/-direction. To do so, we need to 
impose a rotation of —ir/2 rad around the rr-axis. (b) At s = L, we impose two rotations. We 
orient the N th element in the y-direction as in (a) and superimpose a rotation <t>, the control 
parameter, around the y-axis. 



is naturally pointing in the z-direction so that d% is parallel to e 2 , has to be re- 
oriented along the y-direction [see Fig. 7. a]. Mathematically, this translates into 
two functional constraints depending on the positional and rotational degrees 
of freedom, r 1 and q\ alone, 

C° r [r 1 ] = r 1 = 0, (74a) 

T 

= 0. (74b) 

Whereas the first condition (74a) can be physically interpreted as the fixed 
boundary conditions at s = 0, r 1 = 0, the second Eq. (74b) is more diffi- 
cult to interpret due to the lack of direct physical significance of the quater- 
nions. To determine the constraint set in Eq. (74b), we used the relation 
between Euler's principal geometric quantities and a set of unit quaternions, 
whose expression was given in Eq. (10). According to Euler's rotation theorem 
[35], the imposed rotation applied to the first segment (shown schematically in 
Fig. 7(a) can be represented by a rotation angle of — ir/2 rad around the unit 
length vector b = [100] . Following the conversion formula of Eq. (10), the 
equivalent representation in terms of unit quaternions is given by the rotation 
q 1 = [-sin(7r/4) cos(tt/4)] T provided in Eq. (74b). 

At the other end of the rod (s = L), where the rotation is being imposed, 
the last node is fixed at y = d = 22 cm. Furthermore, the last element has 
to be rotated by —pi/ 2 rad around the a;- axis to re-orient along the y- 
direction, as explained above, but we also need to superimpose a rotation $ 
with respect to the y-axis, to simulate writhing [see Fig. 7.(b)]. Mathematically, 
these conditions translate into two constraints depending on the positional and 
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rotational degrees of freedom alone, r +1 and qjv, 
C L r [r N+l ] = r N+1 - [0 d Of = (75a) 
C L q [q N ] =q N - ^[-oos(*/2) sin(*/2) sin($/2) cos($/2)] T = 0. (75b) 

Again, Eq. (75a) is a direct translation of the positional boundary conditions of 
our writhing experiment. Eq. (75b), however, is less intuitive; it is the quater- 
nion representation of the composition of the two rotations shown in Fig. 7.(b). 
The first rotation, A, has already been treated above and can be described by the 
quaternion vector = [— sin(7r/4) cos(7r/4)] T . The second re-orientation, 
B, can be represented by a rotation angle of $ around the unit length vector 
b = [010] T and translates as q^ = [0 sin(<£>/2) cos(<I>/2)] T in the quaternion 
basis. The total rotation imposed to the iV th element is therefore the composi- 
tion of the two rotations A, then B. In terms of quaternions, this rotation q N is 
represented by the multiplication of the two sets of quaternions q 1 ^ and q^ and 
reads q N = q^.q^ which, following the multiplication rule in the quaternion 
basis [35], is given in Eq. (75b). 

To properly account for the boundary conditions in our numerical model 
introduced in Section 4, in addition to the constraints set by Eqs. (74)-(75), we 
also need to include the quantities /x°VC°, ^VCj", A*°VC° and ^VCj, into 
the quadratic vector of equilibrium equations f(u) given in Appendix A. Doing 
so involved the introduction of the Lagrange multipliers fj,®,, fj,^ , fi® and fx q in 
the vector of unknowns u of Eq. (53). Physically, and are the forces 
at each end, written in the Cartesian basis, required to impose the positional 
boundary conditions Eqs. (74a)-(75a) at equilibrium. Similarly and fi q are 
a set of quaternions representing the moments at each extremity, necessary to 
impose the rotational boundary conditions. Finally, these boundary conditions 
can be accounted for in the stability analysis by including the new Lagrange 
parameters and their associated constraints in the second-order condition given 
in Section 4.3 through Eq. (72). 

Finally, before we are able to solve our modified nonlinear algebraic problem 
f(u) = 0, one last important step is needed. The condition given in Eq. (75b) 
is not quadratic in terms of the control parameter and consequently, without 
further modification, the updated nonlinear vector valued function /(it) would 
not be adequate to the numerical framework posed in Section 4. Fortunately, 
there is an appropriate way in the ANM framework to quadratically recast 
Eq. (75b). The technique involves adding two new variables, 

c{a) = cos ($ (a) /2) (76a) 
s(a) = sin ($ (a) /2) , (76b) 

into the vector of unknowns u(a). The procedure is based on the introduction 
of differential equations in terms of the path-parameter a in f(u(a)). Differen- 
tiating Eq. (76) with respect to a, the unknowns (q N , c, s) are now solutions of 
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the quadratic algebraic system, 







q N (a) - — [-c(a) s(a) s(a) c(a)] T , 



(77a) 











dc(a) + -s(a)d$(a), 
ds(a) - ^c(a)d$(a). 



(77b) 



(77c) 



To fully integrate Eqs. (77) into the asymptotic numerical framework of Section 
4.2, we now need to perform a minor modification to the identification technique 
of the power series explained in Eqs. (54)-(58). We recall that the fundamental 
idea behind the ANM is to express the vector of unknowns in a power series of 
a such that, u(a) = it + J2T=i aPu p> whose differential version reads, 



Substituting Eq. (78) into Eq. (77) and identifying the power of a allow us to 
compute the contributions u p of the semi-analytical vector of unknowns u(a), 
following the same procedure described in in Section 4.2. We highlight the fact 
that this method of introducing differential equations in the ANM method is 
a convenient way to represent complex non-polynomial energy functions in our 
numerical model, e.g. to represent highly nonlinear phenomena such as contact 
forces [44]. 

The kinematic boundary conditions, Eqs. (74)-(75), and rotational control 
parameter $, Eqs. (76)-(55), for the writhing problem are all now correctly 
implemented into our algebraic equilibrium equations f(u(a)) = 0, where f(u) 
is given in Eq. (A.l). We can now compute the vector of unknowns, u(a), as 
an asymptotic expansion in terms of the path-parameter, a, using the method 
explained in Section 4, to perform the continuation of the solutions u(a) and 
assess their associated local stability. 

5.3. Comparison between numerics and experiments 

Having introduced, developed and described our theoretical and computa- 
tional tools, we proceed by performing a direct comparison between numerics 
and experiments. In particular, we focus on quantifying the evolution of the 
equilibrium configurations and associated buckling instabilities, as a function of 
the control parameter, $. We highlight that in this comparison, there are no 
fitting parameters; all material and geometric parameters of the experiments are 
independently measured and considered as input variables into the numerics. 

In Figs. 8, 10 and 12b), we compare the top view of some representative 
experimental and numerical equilibrium shapes for a straight rod (ki(s) = 
m -1 ) and a naturally curved rod (ki(s) = 44.84 m _1 ). From these images, we 
measure the maximum transverse displacement of the rod, X max , in the (x,y) 
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(78) 
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Figure 8: Top view of various equilibrium configurations of a straight elastic rod (ki = 
m _1 ) for increasing values of the rotation angle <E>. The experimental pictures have a black 
background and the simulations have a white background. The simulation results are rendered 
to visualize twist by using bi-color rods, (a) Planar equilibrium shape at rest for <E> = 0°. (b) 
Out-of-plane configuration for $ = 720°. (c) Out-of-plane configuration for <J> = 1470°. (d) 
Onset of formation of a plectoneme at the middle of the rod. 



plane (top view) as a function of the rotation angle, $, which is treated as a con- 
trol parameter. Experimentally, the quantity X max was measured from image 
analysis of the digital images taken by the camera located above the apparatus. 
Using these quantities, we then construct the bifurcation diagrams presented in 
Figs. 9(a), 11(a) and 12(a) (for experiments and numerics), for the straight and 
curved rods, respectively. We also analyze the stability of the equilibrium state 
by calculating the first eigenvalue of the stability problem as a function of the 
control parameter, and the results are plotted in Figs. 9(a) and 11(a) (for 
numerics) . In order to quantitatively validate our ANM continuation technique, 
the semi-analytical numerical curves (lines) are superposed onto the experimen- 
tal results (data points, every $ = 30°), for the same value of the control param- 
eter. We highlight, once again, that there are no fitting parameters involved in 
this comparison; all quantities are measured in the experiments, independently 
from the numerics. The excellent quantitative agreement between experiments 
and numerics illustrates the sticking predictive power of our framework. 

We now comment on the experimental and numerical results in more de- 
tail, focusing first on the case of the straight rod (/ti = m _1 ), the results of 
which we plotted in Figs. 8 and 9. Initially, for $ = 0°, the rod exhibits a planar 
equilibrium shape lying in the (y, z) plane due to the effect of gravity. This equi- 
librium configuration is calculated using a classic Newton-Raphson algorithm 
[15] and taken as the initial fixed point ito, the solution of our equilibrium equa- 
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Figure 9: Bifurcation diagram of the straight elastic rod [k\ = m — 1 ). (a) Evolution of 
the maximum transverse displacement, X max = max(abs(X)) (defined in Fig. 8(a)), with 
the control parameter $. Comparison between experimental results and the semi-analytical 
branches computed with MANlab. Solid/dashed lines represent stable/unstable branches, re- 
spectively, (b) Evolution of the first eigenvalue of the stability problem with control parameter 
<I>. The critical angle corresponding to the emergence of plectoneme is $5 im = 1925° for the 
simulations and &% xp = 2025 ± 15° for the experiments. 



tions f(u(a)) = [see Fig. 8. (a)]. When the rotation angle <!> is increased, this 
initial planar shape evolves smoothly into an out-of-plane configuration, sym- 
metric to the (y, z) plane, with an amplitude that grows due to an increasing 
internal twist [see Fig. 8.(b)-(c)]. At a critical value of the rotation angle, <& c , 
the out-of-plane shape loses stability and the rod buckles into a plectoneme state 
[18, 54]: a highly localized structure corresponding to a two-start right-handed 
helix with terminal loops. Beyond this point, our numerical model is no longer 
able to reproduce the rod's configurations since they involve self-contact which 
is not included in our description. To further quantify this process, in Fig. 9 we 
plot the maximum transverse displacement of the rod, X max , as a function of 
the imposed rotation angle, $. Across the full range of $ explored, the exper- 
imental data plotted in Fig. 9 is in excellent quantitative agreement with the 
numerical prediction. 

It is remarkable that the instability threshold for the formation of the plec- 
toneme, $ c , is also well recovered by our local stability analysis showed in 
Fig. 9.(b), where we plot the evolution of the first eigenvalue, <7 1; as a function 
of $. With no fitting parameters, the predicted critical threshold <& c sim = 1925° 
is in excellent agreement (within 5%) with the experimental results <& c exp — 
2025 ± 15°. To plot the semi-analytical bifurcation diagram of Fig. 9, we com- 
puted 80 solution vectors, u(a), expressed in terms of power series expansions 
at the order m = 20, as given in Eq. (55). Using a desktop computer with a 
standard processor (at the time of writing) of 2.71 Ghz and 2.75 Gb of RAM, 
the computation required 9 seconds to determine one asymptotic series. This 
represents a total running time of approximately 12 minutes to simulate the 
full problem, using MANlab. Note that this computational time is mostly due 
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Figure 10: Top view of various equilibrium configurations of a curvy elastic rod (/t i (s) = 44.84 
m — 1 ) for increasing values of the rotating angle, <t>. The experimental pictures have a black 
background and the simulations have a white background. The simulation results are rendered 
to visualize twist by using bi-color rods, (a) Asymmetric out-of-plane equilibrium shape for 
<I> = 150° . (b) One-twist-per-wave configuration for <t> = 900° with three wavelengths between 
the two clamps, (c) One-twist-per-wave configuration for $ = 1800° with five wavelengths 
between the two clamps, (d) Onset of formation of a plectoneme state at one extremity of 
the rod. 

to the 80 inversions of the Jacobian matrix needed to solve the linear systems 
given in Eq. (58), which are of size 1415 x 1415 for N = 100 elements. These 
computations could be made even more efficient by using a dedicated solver such 
as the ones offered by traditional finite element codes but the time optimization 
of the ANM algorithm, which has been investigated [23], is beyond the scope of 
this paper. 

Interestingly, for the case of the naturally curved rod (£i(s) = 44.84 m -1 ) 
the evolution of equilibrium configurations with the rotation angle $ is qualita- 
tively different from the case of the naturally straight elastic rod, as shown in 
Figs. 10, 11 and 12. Once again, the qualitative and quantitative agreement 
between the experimental and numerical results is remarkable. By introducing 
the natural curvature £i(s), the previously symmetric out-of-plane solutions ob- 
tained for the case of straight rods, become asymmetric with respect to the (y, z) 
plane. For small rotation angles the initial planar shape exhibits an asym- 
metric out-of-plane configuration due to the competition between the imposed 
internal twist and the intrinsic twist naturally imposed by k\{s) [see Fig. 10. (a)]. 
Above $ Ri 400°, our results confirm that the rod, jumps into a one-twist-per- 
wave mode due to the presence of natural curvature, as originally reported in 
[55]. In this configuration, the number of waves is equal to the number of twists 
stored in the rod as shown in Figs. 10.(b)-(c). For a critical rotation angle $ c , a 
plectoneme forms, superimposed onto the one-twist-per-wave equilibrium state 
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which is no longer stable. It is interesting to note that, for the naturally curved 
rod, the plectoneme is located at one extremity of the rod rather than at its 
center [see Fig 10. (d)], as found above for the straight rod. 

For the naturally curved rod, our continuation method is able to robustly 
and efficiently follow the equilibrium branches across the full range of consid- 




Figure 11: Bifurcation diagram of the curved elastic rod (k\(s) = 38 m _1 ). (a) Evolution 
of the maximum transverse displacement with the control parameter <t>. Comparison between 
experimental results and the semi-analytical branches computed with MANlab. The presence 
of a constant natural curvature introduces a new wavy configuration with one-twist-per-wave. 
Solid/dashed lines represent stable/unstable branches, respectively, (b) Evolution of the first 
eigenvalue of the stability problem with control parameter ■!>. The critical angle corresponding 
to the emergence of plectoneme is Q c sim = 2919° for the simulations and &% xp = 2895 i 15° 
for the experiments. 
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Figure 12: Hysteretic behavior of the curved elastic rod (rei(s) = 44.84 m — 1 ). (a) Evolution 
of the maximum transverse displacement with the control parameter <t>. Comparison between 
experimental results obtained for increasing and decreasing 4> and the semi-analytical branches 
computed with MANlab. Solid/dashed lines represent stable/unstable branches, respectively, 
(b) The equilibrium configuration for <E> = 150° (two helices with opposite handedness) is 
different from the configuration in Fig. 10(a) for the same control parameter <t>, which is 
significative of multi-stability. 



ered rotation angles, exhibiting excellent agreement with the experimental 
data, as shown in Fig. 11. For this second test case with ki(s) — 44.84 m _1 , 
we used N = 200 elements and we computed 115 asymptotic expansions u(a) 
to determine the full semi-analytical bifurcation diagram, leading to a compu- 
tational time of approximately 35 minutes. Again, the instability threshold for 
the onset of a plectoneme is well recovered by our local stability analysis showed 
in Fig. 11(b). In this case of a more complicated bifurcation diagram, it is re- 
markable that the predicted & c sim = 2919° is within 1% of the experimentally 
measured value of &% xp = 2895 ±15°. Counterintuitively, we find that imparting 
a constant natural curvature to our rods (essentially adding a finite imperfection 
to the stress-free configuration) results in postponing, by approximately 43%, 
the emergence of the plectoneme instability (often synonymous with failure in 
practical systems). To our knowledge, this interesting novel phenomenon has 
thus far been overlooked in the literature and deserves further investigation. A 
systematic study to quantify and rationalize the influence of natural curvature 
on the writhing of a slender elastic rod is beyond the scope of this paper, but is 
an aspect which we plan investigate in future work. 

Another interesting feature in the writhing of the naturally curved rod is 
the hysteric behavior observed for small values of the rotation angle $, before 
entering the one-twist-per-wave regime, as illustrated in Fig. 12. Initially, when 
we increased $, the material twist originally stored due to the natural curva- 
ture Ki(s) is released until the rod jumps to the one-twist-per-wave mode at 
$ w 400°. If we then decrease the observed equilibrium configurations are 
not the asymmetric out-of-plane shapes we previously encountered as illustrated 
in Fig. 10(a). Instead, Fig. 12(b) shows an inversion of helix handedness, known 
as perversion [56] and described as two helices with opposite handedness. If we 
decrease $ even further to negative values, the rod jumps back on the previous 
stable equilibrium branches where the rod twists with the same handedness. In 
the regime of small rotations ($ < 400°), our system is metastable; for the same 
control parameter $, the rod can exhibit two different configurations depending 
on the loading path. Once again, our continuation method correctly predicts 
the different equilibrium states and stability threshold of this complex hysteric 
behavior as shown in Fig. 12. Moreover, this highly nonlinear feature emphasize 
the ability of our numerical technique to follow the equilibrium branches and 
stability of slender elastic rods independently of the complexity of the bifurca- 
tion diagram. 

6. Conclusions and perspectives 

We have presented an original theoretical and computational framework to 
follow the equilibria and stability of slender elastic rods. In our model, we ac- 
count for the elastic energy due to changes of material curvatures and twist, as 
well as the work of external forces and moments, under the assumption that the 
rod is inextensible and unshearable. The main novel feature in our continuation 
method is the use of quaternions to represent rotations. This formulation allows 
for the 3D kinematics to be treated in a geometrically-exact way and result in 
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equilibrium equations that are, at most, quadratic with respect to the state 
variables. We have shown that this quadratic recast of geometric nonlineari- 
ties is particularly well suited for implementation into an asymptotic numerical 
method. This powerful perturbation technique provides access to branches of 
equilibrium solutions in the form of successive portions of power series expan- 
sion by consecutively solving a set of linear systems. The equilibrium branches 
can thereby be followed and their stability evaluated as a function of the con- 
trol parameters. Finally, we have challenged and validated our computational 
framework by considering the specific problem of writhing of a thin rod and 
contrasting our numerical implementation with our own experimental results, 
finding excellent quantitative agreement between the two. We were able to 
successfully and accurately calculate the geometrically-nonlinear configurations 
of the rods, as well as the critical thresholds for instability, with remarkable 
predictive power. We note that our continuation algorithm is able to address 
regions of multi-stability and hysteresis, as in the regime of low rotation angles, 
when increasing or decreasing the control parameter. 

A potential extension of this work would be to incorporate other additional 
mechanical ingredients into our model such as internal stretching, hydrostatic 
loading and contact forces arising, for instance, due to self-contact or when the 
rod interacts with external boundaries. One technical requirement in order to be 
able to introduce new energy terms in our description is that the resulting equi- 
librium equations have to be quadratic to match the present ANM framework. 
However, we have shown that even some some cases of non-polynomial functions 
can easily be reduced to a quadratic form by introducing a limited number of 
new variables in the vector of unknowns in a process we call recasting. Oth- 
erwise, the introduction of these new features can be readily accomplished, as 
long as they derive from a potential energy since continuation methods only 
apply for conservative systems, where an equilibrium can be found. 

We have developed a predictive computational framework to tackle the sim- 
ulations of extreme displacements and rotations in slender elastic rods. Our 
novel method is relatively simple to implement, robust, accurate, flexible and 
computationally efficient. We hope that this technique will be invaluable in 
problems that demand the predictive understanding of the stability, buckling, 
snap-through and other complex mechanical phenomena intrinsic to the extreme 
deformation of slender elastic rods, whose timely revival is highly relevant in a 
variety of currently open problems in both nature and technology. 
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Appendix A. Quadratic form of the vector of equilibrium equations 

The quadratic form of the 14iV-dimensional nonlinear vector valued function 
/ (it) given in Eq. (54) and representing the equilibrium equations (48)-(50) of 
our inextensible and unshearable slender elastic rod, can be written as 
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where the expression of the director d 3 3 (q 2 ), the projection matrix D(q : >) and 
the unknown variables u can be found in Eq. (38), Eq. (47) and Eq. (53), 
respectively. Note that these algebraic equations are the general form of the 
equilibrium equations of the inextensible elastic rod under external forces and 
moments and do not account for particular kinematic boundary conditions or 
control parameter which can vary depending on the problem under study. An 
example of he additional equations needed to impose some fixed boundary con- 
ditions and a rotation angle as control parameter for the case of writhing of a 
rod is given in Section 5. 
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